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We derive exact analytic expressions for the distributions of eigenvalues and singular values for 
the product of an arbitrary number of independent rectangular Gaussian random matrices in the 
limit of large matrix dimensions. We show that they both have power-law behavior at zero and 
determine the corresponding powers. We also propose a heuristic form of finite size corrections to 
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these expressions which very well approximates the distributions for matrices of finite dimensions. 

PACS numbers: 02.50.Cw (Probability theory), 02.70.Uu (Applications of Monte Carlo methods), 05.40.Ca 
(Noise) 
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I. INTRODUCTION 



Spectral analysis of the products of random matrices is a powerful tool in several domains of statistical physics, 
^ ■ allowing, for example, to study Lyapunov exponents for disordered and chaotic dynamical systems |l|. It is also 
useful in a class of problems related to multiplicative matrix-valued noncommutative diffusion processes [2| . Several 
applications go beyond physics, as for instance, those related to the stability analysis of ecological systems J3| or 
to telecommunication applications based on the scattering of electromagnetic waves on random obstacles [J, Q • I n 
many of those cases, some exact analytic results were obtained for relatively small matrices. Interestingly, quite 
often analytic calculations are possible under another limit — the limit of matrix dimensions tending to infinity. 
• • , Examples include products of pseudounitary matrices, representing transfer matrices in mesoscopic wires [6], large 
N Wilson loops in Yang-Mills theory [7H9J or multiplicative diffusion of infinitely large complex and/or Hermitian 
matrices [Ty, [Tl|. In most of these cases, the reason why the exact spectral distribution is within the reach of analytic 
methods is due to a link to free random variable calculus [12lll3|. which is a very powerful technique. This is also why 
the spectra of products of large random matrices represent a challenge for mathematicians [14 , llEfl. In this paper, we 
generalize the analysis of the product of large, square, random Gaussian matrices, performed in 16], to the product 
of rectangular matrices. In particular, we study the product 

P = A!A 2 ...A L (1) 

of L > 1 independent, rectangular, large, random Gaussian matrices A;, I — 1, 2, . . . , L, of dimensions Ni x Ni+i. We 
are interested in the eigenvalue and singular value density of P in the limit Nl + \ — > oo and 

N, 
Ri = —!— = finite, for I = 1,2, . . . ,L + 1. (2) 

Nl+i 
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In other words, all matrix dimensions grow to infinity at fixed rates and, obviously, Rl+i = 1. The product P is a 
matrix of dimensions N% x Nl + \ and has eigenvalues only if it is a square matrix: N\ = Nl+\. 

We assume the matrices A; in the product ([T]) to be complex Gaussian matrices drawn randomly from the ensemble 
defined by the probability measure 

VJV '"' + 1 Tr('A t A, s | 

d/«(A,)«e ^^ l ' j DA,, (3) 

where DA = J\ a 6 d(Re[A] a f,)d(Im[A] a b) is a flat measure. A normalization constant, fixed by the condition 
Jd/x(A) = 1, is omitted. This is the simplest generalization of the Girko-Ginibre ensemble [18l - l20J to rectangu- 
lar matrices. The 07 parameters set the scale for the Gaussian fluctuations in A;'s. The entries of each matrix A/ can 
be viewed as independent centered Gaussian random variables, the variance of the real and imaginary parts being 
proportional to af and inversely proportional to the square root of the number N1N1+1 of elements in the matrix. 

The eigenvalue density of the produ ct liTI) of square Gaussian matrices was calculated in [16| while the singular value 
distribution was determined in 114). Ilq. Il7ll . The eigenvalue density was derived using a planar diagrammatic method 



for non-Hermitian matrices [10|, l2lM23| , while the singular value density was obtained using Free Random Variables 
calculus [13, [24], [25[ . Both techniques work in the infinite matrix size limit. After explaining notation (Section II) and 
listing the main results of the paper (Section III), we shall follow those same methods to derive the corresponding 
results for the product of rectangular matrices. In Section IV, we present a diagrammatic derivation of the moment 
generating function for the product P. In Section V, using the tools of Free Random Variables calculus, we obtain 
the moment generating function for Q = P^P, recovering results given in |17| . Section VI concludes the paper with 
a discussion on particular applications of our results and possible generalizations. 



II. GENERALITIES 

Let us spend a few words on the notations to be used in this paper. The eigenvalue density px(A) of a Hermitian 
matrix X is a real function of real argument, while in the case of a non-Hermitian matrix it is a real function of 
complex argument. In the latter case we shall write px(A, A) and treat A and its conjugate A as independent variables. 

In the Hermitian case, the eigenvalue density can be computed from a Green's function Gx(z) [2y, [231 which 
contains the same information as the density itself: 

Px(A) =-- lim ImG x (A + ie). (4) 

IT e->0+ 

For a non-Hermitian matrix, the corresponding Green's function Gx(z, z) is non-holomorphic and therefore we 
shall write it explicitly as a function of z and z. In this case the eigenvalue distribution is reconstructed from the 
Green's function as [H||j(J 

1 d 

Px{z,z) = - — Gx(z,z). (5) 

7T OZ 

Actually, this equation reduces to (j4j when the non-holomorphic region shrinks to a cut along the real axis, as it 
happens for Hermitian matrices. The Green's function GySz) for a Hermitian matrix is written as a function of a 
single argument since everywhere except on the cut one has 9 ? Gx = 0, and thus it is z-independent. 

In many applications it is often convenient to use the moment generating function, or Af -transform, which is 
closely related to the Green's function: M-^(z) — zGx(-z) — 1. For a Hermitian matrix X one has 

^x(.) = E5 = Ei/^(A)A«dA, (6) 

n>l n>l J 

where the m n 's are the moments of the eigenvalue density. If the matrix X is of finite dimensions N x N , the moments 
are given by m n = -i(TrX n ). The moment generating function encodes the same information as the Green's function 
Gx(z) — z^ 1 Mx.{z) + z^ 1 . Thus, one can calculate the corresponding eigenvalue distribution from Mx(z). 

One can also introduce a similar function for non-Hermitian matrices: Mx_(z,z) = zGx_(z,z) — 1. In this case, 
however, it does not play the role of a moment generating function anymore, since now one also has mixed moments 
(TrX™(X^) fc ), which in general depend on the ordering of X and X^ in the product under the trace. 



The situation is slightly simplified when the M-transform is a spherically symmetric function: Mx(z,z) — 
■Mx(\z\ 2 )- In this case equation §5§ can be cast into the form 

p^z,z)^-M' x (\z\ 2 ) + f5 2 (z,z) (7) 

where VW X is the first derivative of .Mx and / = 1 + A^x(O) is a constant representing the fraction of zero modes. In 
this case, the eigenvalue distribution is spherically symmetric as well (see for example [16] for the product of square 
matrices). As we shall see later, this is also going to be the case for the product ([l} of rectangular Gaussian matrices 



III. RESULTS 

The matrix P ([T]) has eigenvalues only if it is square, while it has singular values for any rectangular shape. As a 
matter of fact, its singular values can be determined as the square roots of the non-zero eigenvalues of the matrix 

Q = P f P (8) 

or, alternatively, of the matrix R = PP^. Q and R are Hermitian, and they have non-negative spectra which differ 
only in the zero modes. 

The main finding of this paper is that the eigenvalue distribution and the M-transform of the product |T]) are 
spherically symmetric. We shall show the M-transform to satisfy the L-th order polynomial equation 

where the scale parameter is a = o\<J2 ■ ■ ■ &l- When all of the matrices involved are square, this equation reproduces 
the results in flol. 



An analogous equation for Q reads 



f - M Q (z) + l A( M Q (z) \ z 

^T^llt-RT + 1 J-^ ( 10 ) 

It is an algebraic equation of order (L+l), and it was first obtained in [17| in the context of wireless telecommunication. 
Equations © and (fit)]) are strikingly similar. They actually differ only by the prefactor in front of the product. 
Moreover, the free argument in the first equation is |z| 2 , while z in the second one. This observation represents the 
second main result of this paper. Since conventions used in telecommunication theory and in physics differ a bit, in 
Section V we rederive equation (JTUJ) for completeness. 

When P is a square matrix, then Rx = 1 and the square root at the beginning of equation (fTU|) can be omitted. 
When the product of square matrices is considered, all of the i?/'s become equal to unity and the two equations take 
the following form: 

(AMM 2 ) + 1) L = ^ , M Q \z)(M Q (z) + l) L+1 = ^. (11) 

Equations (jH]) can be easily rewritten in terms of the corresponding Green's functions (see the previous section). If 
one does that and then applies the prescriptions in ([3]) and (U) respectively, it becomes clear that 

pp(A,A) - lAI" 2 ^, and p Q (\) ~ \~ ttt : as A -> 0. (12) 



In the more general case of rectangular matrices, when solving equations (|9]) and (|T0|) for the Green's functions, one 
can then see that only those brackets in which Ri — 1 contribute to the singularity at zero, while all others approach 
a constant for z —¥ 0. Thus, the eigenvalue density displays the following singularity 

pp(A,A)~ |Ar 2£ ^\ as A^O, (13) 



where s is the number of those ratios among R± , 
eigenvalue density of Q behaves as 



. , Rl which are exactly equal to unity. On the other hand, the 



/Oq(A)~A .+i, 



A^O. 



(14) 



The third result we want to mention here is a heuristic form for the finite size corrections to the eigenvalue 
distribution. For a large but finite order of magnitude N of the matrices involved, the eigenvalue distribution is still 
spherically symmetric. So, let pn(j~) denote the radial profile of this distribution, where r = |A|. As we shall show, 
the evolution of the radial shape with the size N is very well described by a simple multiplicative correction: 



Pn{t) = p{r)-eiic (q(r - ct)Vn) 



(15) 



In the N — > oo limit the correction becomes a step function, so that Poo{f) = p{r) for r < a and Poo(r) — otherwise. 
The shape of the limiting radial distribution p(r) comes from the solution of (|S]). This type of finite size corrections 
can be derived analytically for Girko-Ginibre matrices |3ll - l33| . Here we show that it also works very well for the 
eigenvalues of the product of Gaussian matrices. It is very generic and possibly applies to other random matrix 
ensembles with spherically symmetric eigenvalue densities. 



IV. THE EIGENVALUES OF A PRODUCT OF RECTANGULAR GAUSSIAN RANDOM MATRICES 

In this section, we present a derivation of the main result of this article, equation (j9|), a realization of it in the 
case L = 2, and numerical simulations to confirm our findings. To this end, we employ a technique for summing 
planar diagrams, the Dyson-Schwinger equations (more precisely described in [Hi, LL6|), extended to a non-Hermitian 
framework. 

The evaluation of the Green's function for a product P of random ensembles by means of planar diagrammatics 
is non-linear w.r.t. the constituent matrices. It is possible to linearize the problem by means of the following 
trick [nj, LL6[ ■ Consider the following block matrix: 
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It is a matrix of dimensions N to t. X -Aftot.j where N tot . 
power of P is a block diagonal matrix [ljj, [l6| 









Ni + N2 + . . . + Nl ■ It is important to notice that the L-th 
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(17) 



bJ 



with square blocks Bi = A1A2 . . . Ax,_iAl, B2 = A2...A^Ai, ... , being cyclicly permuted products of Ai, 
A2, . . . ,Al. All these blocks have identical non-zero eigenvalues. They differ only in the number of eigenvalues, 
which may vary from block to block. The first diagonal block Bi is equal to the product P ((TJ). Taking into account 
that this block has dimensions NixNi, while the total matrix P L has dimensions AT to t. X AT to t., one can easily deduce 
the following relation between the M-transforms of P and P: 



M p (w,w) = -^tf P (w L ,w L ) 

-IVtnt. 



(18) 



The importance of this relation relies in the fact that one can use it to calculate Mp(z, z) from Mp(w, w). The latter 
can be calculated using Dyson-Schwinger equations, since the matrix P is linear w.r.t. the constituent matrices Aj. 



The first step-in writing the Dyson-Schwinger equations is to know the propagators of the random matrix in 
question, namely P, or more precisely its "duplicated" version: 



P 
P t 



(19) 



We shall think of it as a four-block matrix, each block being anixi block matrix. We shall denote the L x L 
block indices in these four blocks by Im (upper left cornerV Im (upper right), Im (lower left), Im (lower right), each 
one covering the range 1,2, . . . , L; for example [P D ]2i = AJ\ All the other matrices involved shall inherit this same 
structure. For instance, 



G D = 



G 



G< 



G< 



G< 



n cD ]n [ G °]i 2 

[ G °] 21 [ G °] 22 

[ gD 1m ^ d U 



[ GD ]n [ GD ]t 2 
[G D ] 21 [G°] 22 



V [G»]- L1 [G°] 



L2 



[ GD L 



J2i 



m^ 



[G D ' 



IlL 



2L 



[G D ] 



LL 



[G D ]n [G°] 12 
[ G °] 2 T [ GD ] 2 , 



G 1 



ML. 



[G D 1 



[G D ]tt [&]- 2 
[G D ]n [G D ]^ 



[G D ] lX \ 
[G D 



G 1 



J2L 



1 T.T. 



[G D Jtt 
[G°k 



[G D ]rr/ 



(20) 



and similarly for W D and S D (to be defined in a moment). For the sake of simplicity, we shall disregard some 
subscripts and symbols of dependence on w, w. 



We are interested in computing the Green's function of P, i.e. 

L 

i i 

Gp(w,w) 



rp r Q WW . 

Ntot. iVtot 



tot. l=1 iVtot. l=l 



NiQu, 



where it is useful to define the normalized traces 



0* - ^ [G°], . Qu - ^Tr [G»] a , Q- u . i-Tr [G"] B , Q Tl - ±-Tr [G°] n 

Hence, we should evaluate the Qua. 

The only non-zero propagators of P D are readily determined from the probability measures in (|3]): 
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21/ y/NiN; 
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32/ y/TWh' 



=ljVi <E>1jvi- 



(23) 



11 - JlX/ y/N L Ni 

Thus, we are now in position to write down the two Dyson-Schwinger equations for P D . The first one, being the 



definition of the self-energy matrix S D , is independent of the propagators 

G D = (W D -S D ) _1 



squat 

m 



(24) 



where W D is defined as ioljv tot in its left upp er block (where w £ C), Wljv tot in the right lower block, and zero 
elsewhere. The second one is presented in [10|, LL6|, and the structure of the propagators (l23l) implies that the only 
non-zero blocks of the self-energy matrix read 



[s d Lt = 



vW 



1+1 



= Tr [G L+ij+r 1 N l =°i 



Ni 



+i 



Ni 



Si+ii+i 1 ^, 



(25) 



[ED) » " 7W^w, T ' [QD] ^'-' lw " ""Vt^- ^ 



W- 



(26) 



for all Z = 1, 2, . . . , L, with the cyclic convention = L, where the normalized traces (|22|) have been used. 

Results Ij25|> . (|26| mean that the four blocks of the matrix (W D — S D ) are diagonal. Such a matrix can be 
straightforwardly inverted: its four blocks remain diagonal, and read 



w7 2 1jv 2 



(W D - S r 
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(27) 



wj l 1 Nl I 



where, for all I = 1, 2, . . . , L, 

1 



i 12 a i i2 / \2 \/Ni-iNi + i 



(28) 



Substituting (f2"T|) into (|2"3|) , we find that the only non-zero blocks of the duplicated Green's function (j2"0)) are, for all 
1 = 1,2,. ..,L, 

[G D ] n =W7/l Wi , [G D ] a = aj7Jl^„ [6°]^ = ^!^, [G D ] n = Wll l Nl . (29) 

Taking the normalized traces of both sides of every equality in f|29[) leads to the final set of equations, 

Gu = wji, Q a = a/7;, Q- u = firfu Gn = W "U- ( 30 ) 

The structure of equations (I3U1) is the following: the fourth one is the conjugate of the first, and it is then 
redundant. The second and third ones read 



c ^ lNl 



Hi - u i 



N, 



-Gi+ij+ilu 



(31) 



2 Nl-l r 

yji = a i-i\j -j^yj^ij-iii- 



(32) 



We see that (j28|) , (131J) and (|32|) form a closed set of 3L equations for 3L unknowns, 5,7, C/y, and 7;. Once solved, when 
the 7/'s have been found, we are able to recover the Green's function of P, and subsequently the M-transforms of P 
and P (in the argument w L ) (|18[) . 



1 L _ 1 L 

G=(w,w) =w — — ^Nai, i.e., Mp(w,w) = — — ^JJVj/ij, i.e., 



1 



""*' J=l """■ Z=l 

where we have traded the 7;'s for a more convenient set of variables, 

^ = |w| 2 7; - 1. 



;=i 



(33) 

(34) 



Equations (j3"Tj) and (13"2l form a set of decoupled recurrence relations for Qq and £/j ; , respectively. Iterating these 
recurrences down to I — 1 gives us: 



Gil = Gn 



1 1 



(criCr 2 ...CTi_i) 2 V / ^7i72---7/-i' 



(35) 



2 1 

Qii = Gii(vi02---m~i) -7=72... n- (36) 

Applying the cyclic convention = L, we get the following equation: 

(cricr 2 ■•■o'l) (7172 • • • 7i) Su = G X J- (37) 

Straightforwardly, we get a trivial solution: Qjj — for all I, i.e., remembering (1281) . 7/ = l/|w| 2 , or equivalently //; = 
from ([34l) . and therefore Mp(z,z) = from ([33]) . This is the holomorphic solution, holding outside the eigenvalue 
density domain. In order to retrieve information on the eigenvalue distribution, let us take Q^j 7^ 0. 



After a change of variables to \ii and some simplifications, (|2"8j) becomes 

Mi + 1 Ri ' 



„ , H%fa ' , (38) 



from which we get the relation: 



Ri Ri 

After plugging (1591 into (|37p , in terms of the /!/ variables we obtain 






"" h| £ + 1 )- £ + 1) — 



On the other hand, substituting (1391) into (1331) . we get 

M-p (w L ,w L ) = m. (41) 

All in all, after changing the argument from w to z = w L we see that Mp(z, z) obeys the L-th order polynomial 
equation 

M^ + 1 )(*^ + 1 )(MZ^ + 1 )J4, (42) 

i?i J \ R 2 J \ R L J a 2 

which is precisely the first main result of our article, ([9]). 

The last point to be addressed is to determine the validity domain of the non-holomorphic solution (|42j) . knowing 
[22| that on the boundary of such a domain, the non-holomorphic and holomorphic solutions must be joined. Thus, 
plugging the latter (Mp(z,z) = 0) into (|42|) . we obtain an equation for the borderline: 

\z\ - o. (43) 

This means that the eigenvalues of the P matrix are scattered on average, with the density stemming from (|42p . 
within a centered circle of radius a. 

When L = 2, (J42"T) is just a second degree equation, and it is easily solved. Indeed, in this case the non-holomorphic 
M-transform reads 



\[-l-R+\J(l-R) 2 + ±R^ 2 



M P (z, z) = - -1 - R+ W(l - R)* + AW-4- , (44) 



where we pose R = R2 = N2/N1, and where the proper solution of (|4"2"j) has been picked up in order to satisfy the 
matching condition (|43[) with the holomorphic one on the borderline. As a result, we immediately obtain the Green's 
function: 



1 1-^(1-^+4^ 



G P (z,z) = — ( 1-R+\I(1-R)* + 4RL-L.\ . (45) 



8 

When deriving the average spectral density, one has to be cautious in the vicinity of the origin of the complex 
plane in order to properly take possible zero modes into account. Let us first expand (|45|) near z — in order to 
clarify its behavior: 

Gp(z,z) ~ — h regular terms, as z — > 0, where /=< ' ' (46) 

z I 0, for R > 1. 

Taking the derivative {l/ir)dz of this singular term yields a Dirac delta function at the origin, f8^ 2 \z, z). Alto- 
gether, 

\ Ik*-, f ^T + f^K^), for \z\<a, 

p P (z,j) = { y/(l-R)>+4Rlg (47) 

0, for \z\ > a. 

Moreover, one can also verify that the density, in the thermodynamic limit, changes on the borderline from being 
non-holomorphic with value given by 

1 1 L 1 

pp(z,z) Li =0 . = — ji?h, where "5" = X! ~5~' ( 48 ) 

■no 1 R h *-^ Ri 

to being holomorphic with value 0. However, for finite sizes of the random matrices, this step gets smoothed out. Let 
us then consider the radial density, 

p¥ d -(r) = 2irrp-p(z,J)\ lzl=r (49) 

and, following lGj, let us propose the following model for this finite- N effect (where by N we denote the order of 
magnitude of the dimensions of the matrices, say N = N\). We introduce an "effective" radial density in order to 
properly incorporate such finite-iV behavior at the borderline, 

pf- (r) = / P ad - (r) ierfc (q(r - a) Vn) , (50) 

where q is a free parameter whose value is to be adjusted by fitting. We numerically verify this hypothesis (see figures 
[T] and I2J). 



V. THE SINGULAR VALUES OF A PRODUCT OF RECTANGULAR GAUSSIAN RANDOM 

MATRICES 

In the following we show how to derive formula (|10p. i.e., an (L + l)-th order polynomial equation obeyed by 
the M-transform (which, as already discussed, encodes the same information contained in the spectral density) of 
the Hermitian matrix Q = ptp (jHJ), P being the product ||T]) of rectangular fl2J) Gaussian random matrices ©• Our 
result agrees with that in |17j , obtained in the context of wireless telecommunication theory, provided we synchronize 
the conventions. In particular, our resolvent G(z) relates to Stieltjes transform as G{z) = —G(—s). The underlying 
idea will be to rewrite Q as a product of some Hermitian matrices in order to apply the techniques provided by Free 
Random Variables (FRV) calculus. Loosely speaking, FRV calculus (initiated by the pioneering works of Speicher 
and Voiculescu et al.) can be thought as the extension of standard probability theory to non-commutative objects, 
such as matrices. Given the broadness of the topic, we shall not attempt any introductory discussion here, and we 
refer the non-expert reader to [lj, [3 • 

Let us commence by defining, for any I = 1, 2, . . . , L, a square -/V; + i x Ni + i matrix 

Q ; = (AiA 2 . . . A^A,) 1 (AiA a . . . A^A,) = AjAJ^ . . . A^A t 1 A 1 A 2 . . . A,_iA,, (51) 

being a generalization of Q which includes only the first I random matrices, as well as a square Ni x Ni matrix, which 
differs from Q; only in the position of the last matrix in the string, i.e., A;, which is now placed as the first matrix 
in the string, 

Q, = A ; A| Aj_ a . . . A| At AiAa . . . A_ x = (a ; A}) Q,.,, (52) 
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FIG. 1: Numerical verification of the theoretical formula (|47[) for (the radial part (|49[) of) the mean spectral density pp(z, z) 
of the product P of L — 2 rectangular Gaussian random matrices, as well as the finite-size correction (|50l) . 
UP LEFT: A numerical histogram (the black line) versus the theoretical prediction (1471) . supplemented with the finite-size 
smoothing ([50]) (the red plot), for Nt = 100 and N 2 = 200 (i.e., R= R 2 = 2), and for 10 5 Monte-Carlo iterations (i.e., the 
histogram is made of 10 eigenvalues). The adjustable parameter q (|50[) is fitted to be q w 1.14. 

UP RIGHT: An analogous graph to UP LEFT, this time with iVi = 100 and N 2 = 150 (i.e., R = 1.5). We find q « 1.08 here. 
DOWN LEFT: An analysis of the finite-size effects: numerical histograms for Ni - 50, N 2 - 100 (black), Ni = 100, N 2 = 200 
(dashed red), N\ — 200, 7^2 = 400 (dotted blue), i.e., with the same rectangularity ratio R = 2, but increasing matrix dimen- 
sions. We observe how these plots approach the green line of the theoretical formula (|47[) for the density in the thermodynamic 
limit. 

DOWN RIGHT: Numerical histograms for the matrix sizes of Ni = 100, N 2 = 200 (i.e., R = 2; black) and Ni = 200, N 2 = 100 
(i.e., R — 1/2; red). Due to the presence of the zero modes (not displayed in the picture), the latter is half of the former. 



We are interested in the eigenvalues of the Hermitian matrix Q = Q^ . 



The orders of the terms in the two above products (|5T|) . (1521) are related to each other by a cyclic shift, therefore, 
for any integer n > 1, there will be TrQ™ = TrQ". Hence, the M-transforms (see equation ([5])) of the two above 
random matrices are related by the following relation 



M Ql (z) 



V — (TrQI 1 ) 

f^ z« N l+1 y ^ l ' 



n>\ 



N l+1 ^z"N^ " 

n>l 



Ri 



R, 



i+i 



M & (z). 



(53) 



Now, let us consider the functional inverse of the M-transform, called the iV-transform, defined as: Mq, (JVqj (2)) 
A^q, (A/q, (z)) = z. Employing this definition within equation (|53[) one easily obtains 



«*(•)-"* ^ 



Ri-\ 



(54) 
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FIG. 2: Analogous graphs to figure [Q DOWN LEFT, but for L = 3 (LEFT) and L = 4 (RIGHT). 



Now, since it can be safely stated that independent random matrices become free with respect to each other 
in the thermodynamical limit, it becomes clear that the reason for introducing the auxiliary matrix Q; is that it 
is a product of two free matrices, A/AJ and Q; 1. Then, the FRV multiplication [l3| law for free matrices can be 
applied. Such law states that the A-transform of the product of two free matrices, A and B, is simply given by 
Nab{z) = 2/(1+2) Na(z)N-b(z). (In the language more often found in the literature on the subject, the A^-transform 
is replaced by the so-called S'-transform, Sx(z) = (z + 1) / (z Nx(z)) , which then obeys a simpler multiplication law, 
Sxb(z) = S\(z)S-b(z)). So, when applying this relation to the Q/ matrix (|52[) one can write, for I = 2, 3, . . . ,L 



n q,^ = VTT NA ^ iz)NQ '- l{z) - 



(55) 



From equations (I54[) and (|55[) . we now eliminate the Af-transform of the auxiliary Q;, which leaves us with the 
following recurrence relation for the A-transform of Q;, 



N Ql (z) 
with the initial condition, 



z + 



Ri JV A,Aj 



Ri 



+1 



Ri 



z)N, 



Qi-i 



R 



l+l 



Ri 



for 



( — 2, 3, . . . , L, 



N Ql {z) = N^ 



R% 

Ri 



N ^M (|* 



(56) 



(57) 



which stems from (J54I) and from the fact that Qi = AiAJ. The solution of this recurrence (f56|). (|5T[) is then readily 
found to be 



#Q t (*) 



V L-1 



(z + R 2 ) (z + R 3 )...(z + R L ) N ^ A l [rJ Na * a $ {r 2 J ■ --^At [ R J 



(58) 



,t 



It remains now to find the A^-transforms of the random matrices A; A]. They are examples of the so-called 
"Wishart ensembles", and the problem of computing their A-transforms, with the same normalization of the proba- 
bility measures (O of the A/'s which we are employing, has first been solved in [3J|: expressions (1.8), (2.8), (2.13), 
(2.14) of this article yield the Green's function of AjAl, which immediately leads to the pertinent A^-transform, 



VW = ff ' 



(* + 1) (v^&+ \fW) 



Substituting (|59| into (|58jl. one finally arrives at the desired formula for the A"— transform of Q = Ql, 

I I ■ • . I t- + 1 



N Q (z) = a 2 ^Ih-(z + 1) ( -J" + 1 
z \Hi 



z 

R~2 



\Rl 



(59) 



(60) 
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with a defined as in the previous sections. In other words, the corresponding M-transform Mq(z) satisfies the 
following polynomial equation of order (L + 1), 



[jr 1 (M M , u ( Mq(z) \( M Q (z) \ ( M Q (z) \ z 

v Rl md {Mq{z) + 1} {-or + v {-^r +1 )- {-nr +1 )~^ 

or in the case of Nl+i — N\ (i.e., R\ = 1, required when one wishes for P to have eigenvalues too), 



Mq(z) 
This completes our derivation of ((TU 



' (M Q (,) + W^M + 1 



V #2 



M Q {z 
Rl 



(61) 
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FIG. 3: Numerical verification of the theoreticai formula (|61|l for the mean spectral density pq(A) of the random matrix 

Q = P^P ([H|. Everywhere we have Nl+i — N\ — 50. The number of Monte-Carlo iterations is 20,000, i.e., all the histograms 

are generated from 10 6 eigenvalues. 

LEFT: L = 2, and the matrix sizes are chosen to be N\ = 50, N 2 = 150. 

MIDDLE: L = 3, and the matrix sizes are Ni = 50, N 2 = 100, N 3 = 150. 

RIGHT: L = 4, and the matrix sizes are iVi = 50, N 2 = 100, iV 3 = 150, N 4 = 200. 

We have performed extended numerical tests of the formula (|61[) . in all cases obtaining perfect agreement, see 
figure [31 



VI. CONCLUSIONS 



The main contribution of this article is equation (|9]) for the M-transforms of the product P = Ai A2 . . . Al (H|) 
of an arbitrary number L of independent rectangular ((21) Gaussian random matrices ([3]) . Knowing the M-transform 
one can easily calculate the eigenvalue density of the product (JT]) , which turns out to be spherically symmetric in the 
complex plane. We also discussed a striking resemblance of equation ([9]) to the corresponding equation (p~0|) of the 
Hermitian matrix Q — P^P IJ5J), whose eigenvalues are equal to the squared singular values of P. Both these equations 
are polynomial (of orders L and (L + 1) respectively), so in general they may only be solved numerically; however, 
some properties of the mean spectral densities can still be retrieved analytically, such as their singular behavior at 
zero dI5b. (HI). 



We are tempted to conjecture that this similarity of the M-transforms for P and Q is generic for random matrices 
possessing rotationally symmetric average distribution of the eigenvalues, and that the corresponding equations differ 
only by the prefactor which we have discussed while comparing ([9]) and (| 10[) . For such models, the non-holomorphic 
M-transform Mx.(z,"z) is a function of the real argument \z\ 2 , thereby allowing for functional inversion, and hence for a 
definition of the "rotationally-symmetric non-holomorphic iV-transform" — even though for general non-Hermitian 
random matrices a construction of a "non-holomorphic iV-transform" remains thus far unknown. This new ./V- 
transform is then conjectured to be in a simple relation to the (usual) ./V-transform of the Hermitian ensemble X^X. 
In a typical situation, the latter will be much more easily solvable than the former, owing to the plethora of tools 
devised in the Hermitian world, albeit the opposite may be true as well. This is indeed the case here — our derivation 
of ([5]), based on non-Hermitian planar diagrammatics and Dyson-Schwinger's equations, is much more involved than 
a simple application of the FRV multiplication rule leading to (jlOp — and consequently, the aforementioned hypothesis 
would provide a shortcut to avoid complicated diagrammatics. To the best of our knowledge, this would be the first 
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use of Free Random Variables calculus to compute the mean spectral density of a non-Hermitian product of random 
matrices. 

We have also suggested a heuristic model of the finite-size behavior of the density of P near the edge of the 
eigenvalues support (|5P|) . deducing it from analogous considerations |3ll-l33j made for the Girko-Ginibre ensemble, 
where this behavior is known analytically. It performs outstandingly well when checked against numerical simulations. 

Let us also remark that one could argue, as for square matrices, that the large- N limit result is the same for elliptic 
Gaussian ensembles [16j ■ We also believe that one can further weaken the assumptions on the matrices involved, just 
requiring them to belong to the Gaussian universality class of matrices having independent entries and fulfilling the 
Pastur-Lindeberg condition [38] (the matrix analogue of the generalized central limit theorem in classical probability 
theory [39] ) . One unexpected implication of such universality is that a product of random matrices whose spectra do 
not necessarily display rotational symmetry has an eigenvalue distribution which does possess rotational symmetry 
on the complex plane (i.e., the average density depends only on |A|). 

Let us now list some possible applications of these results to wireless telecommunication, quantum entanglement 
and multivariate statistical analysis. 

Information theory for wireless telecommunication has been intensively developed in the past decade, after it had 
been realized that in a number of situations the information transmission rate can be increased by an introduction of 
multiple antenna channels, known as the "multiple-input, multiple-output" (MIMO) transmission links. The MIMO 
capacity for Gaussian channels has been calculated in the pioneering work [351 ] . triggering large activity in the field. 
Immediately, it became clear that an appropriate language and methods to address this type of problems are provided 
by random matrix theory (consult [5] for a review) . The model considered in our paper can be applied to a situation 
of signals traveling over L consecutive MIMO links. The signal is first sent from Aq transmitters via a MIMO link 
to A^ receivers, which then re transmit it via a new MIMO link to the subsequent N3 receivers, etc. Clearly, the 
capacity will depend on these numbers of intermediate re-transmitters; in particular, if any of the Ni's is small, 
the capacity will be reduced. The effective propagation is given by the matrix P = Al . . . A2A1. Such a model of 
multifold scattering per propagation path has been already proposed in [171 ] , where the moment generating function, 
the M-transform, for P^P was calculated. Our result for the M-transform for P complements this calculation. 

Let us also mention that one could imagine a more general situation, where MIMO links form a directed network 

- each directed link Im representing a single MIMO channel between TV; transmitters and N m receivers. (The 

previously discussed case corresponds to a linear graph, 1 — >• 2 — >• . . . — >• L.) A complex directed network of MIMO 

links is somewhat similar to the structures appearing in the context of quantum entanglement. There, one considers 

graphs whose edges describe bi-partite maximally entangled states, while vertices describe the couplings between 



subsystems residing at the same vertex 36j. In the simplest case of a graph consisting of a single link, it is just a 
bi-partite entangled state. The corresponding density matrix for a bi-partite subsystem is given by Q = A^A, where 
A is a rectangular matrix defining a pure state, being a combination of the basis states in the subsystem, \a a ) and |/3&) 
(see for instance [20]). One can easily find that linear graphs with additional loops at the end vertices correspond. The 
density matrix for the subsystem sitting in the end vertex is given by Q = P^P, where P = A1A2 . . . A^ [36]. If all 
the subsystems are of the same size, the average spectral distributions [1J, [l5| of Q are known as the "Fuss-Catalan 
family" [32}; they can be obtained from PU|) by setting all the i?/'s to 1. However, if the subsystems have different 
sizes, one needs to apply our general formula (flU|) . 

Finally, another area of applications of our approach is related to multivariate analysis. The main building block 
there is the Wishart ensemble, corresponding to L = 1 in our formalism. The link between the spectral properties of 
P and Q may allow one to avoid the well-known bottleneck caused by the non-Hermiticity of time-lagging correlation 
functions. This issue will be discussed in a forthcoming publication. 
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